library(readxl)
mydata <- read_xlsx("Bank.xlsx")
mydata <- as.data.frame(mydata)
head(mydata)
##   Client_id Credit_score Gender Age Tenure Products_number
## 1  15662021          685 Female  42      2               2
## 2  15634359          639 Female  41      5               1
## 3  15697270          608   Male  27      4               1
## 4  15762331          569   Male  37      9               1
## 5  15709136          620   Male  28      8               2
## 6  15741719          540 Female  40      3               1
##   Credit_card Active_client Churn
## 1           0             0     0
## 2           1             0     0
## 3           1             1     0
## 4           1             0     0
## 5           1             1     0
## 6           0             1     0

Description:

mydata$Gender <- as.factor(mydata$Gender)

mydata$Credit_card <- factor(mydata$Credit_card, 
                            levels = c("0", "1"), 
                            labels = c("No", "Yes")) 

mydata$Active <- factor(mydata$Active_client, 
                        levels = c("0", "1"), 
                       labels = c("No", "Yes")) 

mydata$Churn <- factor(mydata$Churn, 
                       levels = c("0", "1"), 
                       labels = c("No", "Yes")) 
summary(mydata[c("Churn", "Credit_score", "Gender", "Age", "Tenure", "Products_number", "Credit_card", "Active")])
##  Churn       Credit_score      Gender          Age       
##  No :7963   Min.   :350.0   Female:4543   Min.   :18.00  
##  Yes:2037   1st Qu.:584.0   Male  :5457   1st Qu.:32.00  
##             Median :652.0                 Median :37.00  
##             Mean   :650.5                 Mean   :38.92  
##             3rd Qu.:718.0                 3rd Qu.:44.00  
##             Max.   :850.0                 Max.   :92.00  
##      Tenure       Products_number Credit_card Active    
##  Min.   : 0.000   Min.   :1.00    No :2945    No :4849  
##  1st Qu.: 3.000   1st Qu.:1.00    Yes:7055    Yes:5151  
##  Median : 5.000   Median :1.00                          
##  Mean   : 5.013   Mean   :1.53                          
##  3rd Qu.: 7.000   3rd Qu.:2.00                          
##  Max.   :10.000   Max.   :4.00
fit0 <- glm(Churn ~ 1, 
            family = binomial, 
            data = mydata)

summary(fit0)
## 
## Call:
## glm(formula = Churn ~ 1, family = binomial, data = mydata)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.36333    0.02483  -54.91   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 10110  on 9999  degrees of freedom
## Residual deviance: 10110  on 9999  degrees of freedom
## AIC: 10112
## 
## Number of Fisher Scoring iterations: 4
fit1 <- glm(Churn ~ Credit_score + Gender + Age + Tenure + Products_number + Credit_card + Active,  
            family = binomial, 
            data = mydata)
library(car)
## Loading required package: carData
vif(fit1)
##    Credit_score          Gender             Age          Tenure 
##        1.000914        1.003033        1.079058        1.001690 
## Products_number     Credit_card          Active 
##        1.001769        1.000882        1.078969
mydata$StdResid <- rstandard(fit1)
hist(mydata$StdResid,
     main = "Histogram of standardised residuals",
     ylab = "Frequency",
     xlab = "Standardised residuals")

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# Calculate Cook's distance
mydata$CooksD <- cooks.distance(fit1)

# Define Cook's distance threshold
cd_threshold <- 4/nrow(mydata)

# Display apartments with Cook's distance above the threshold
influential_units <- mydata %>%
  filter(CooksD > cd_threshold) %>%
  arrange(desc(CooksD)) %>%
  select(Client_id, Churn, Credit_score, Gender, Age, Tenure,
         Products_number, Credit_card, Active, StdResid, CooksD)

head(influential_units, 3) %>%
  mutate(CooksD = round(CooksD, 4))
##   Client_id Churn Credit_score Gender Age Tenure Products_number
## 1  15673591   Yes          842   Male  44      3               4
## 2  15664754   Yes          640   Male  39      9               4
## 3  15805704   Yes          745 Female  32      2               4
##   Credit_card Active StdResid CooksD
## 1          No    Yes 3.492270 0.0022
## 2          No    Yes 4.088059 0.0021
## 3          No    Yes 3.975871 0.0021
anova(fit0, fit1, test = "Chi") 
## Analysis of Deviance Table
## 
## Model 1: Churn ~ 1
## Model 2: Churn ~ Credit_score + Gender + Age + Tenure + Products_number + 
##     Credit_card + Active
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1      9999    10109.8                          
## 2      9992     8831.1  7   1278.7 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
logLik_fit0 <- as.numeric(logLik(fit0))
logLik_fit1 <- as.numeric(logLik(fit1))

McFadden_R2 <- 1 - (logLik_fit1 / logLik_fit0)
McFadden_R2
## [1] 0.1264813
summary(fit1)
## 
## Call:
## glm(formula = Churn ~ Credit_score + Gender + Age + Tenure + 
##     Products_number + Credit_card + Active, family = binomial, 
##     data = mydata)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -2.7841506  0.2273658 -12.245  < 2e-16 ***
## Credit_score    -0.0006461  0.0002756  -2.345 0.019052 *  
## GenderMale      -0.5313055  0.0535313  -9.925  < 2e-16 ***
## Age              0.0723159  0.0025305  28.577  < 2e-16 ***
## Tenure          -0.0140652  0.0092067  -1.528 0.126584    
## Products_number -0.1755235  0.0452403  -3.880 0.000105 ***
## Credit_cardYes  -0.0400537  0.0582808  -0.687 0.491923    
## ActiveYes       -1.0777731  0.0568705 -18.951  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 10109.8  on 9999  degrees of freedom
## Residual deviance:  8831.1  on 9992  degrees of freedom
## AIC: 8847.1
## 
## Number of Fisher Scoring iterations: 5
exp(cbind(OR = fit1$coefficients, confint.default(fit1))) 
##                         OR      2.5 %     97.5 %
## (Intercept)     0.06178154 0.03956622 0.09647013
## Credit_score    0.99935406 0.99881439 0.99989402
## GenderMale      0.58783702 0.52928669 0.65286426
## Age             1.07499485 1.06967635 1.08033979
## Tenure          0.98603323 0.96839999 1.00398754
## Products_number 0.83901769 0.76782543 0.91681085
## Credit_cardYes  0.96073788 0.85703032 1.07699489
## ActiveYes       0.34035262 0.30445343 0.38048482
mydata$EstProb <- fitted(fit1)

mydata$Classification <- ifelse(
  test = mydata$EstProb < 0.50,
  yes = "No",
  no = "Yes"
)

mydata$Classification <- factor(
  mydata$Classification,
  levels = c("No", "Yes")
)

Table <- table(
  Actual = mydata$Churn,
  Predicted = mydata$Classification
)

Table
##       Predicted
## Actual   No  Yes
##    No  7756  207
##    Yes 1709  328
# Install once if not already installed
# install.packages("pROC")

library(pROC)
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
# Loads the pROC package, which provides functions to compute and plot
# Receiver Operating Characteristic (ROC) curves and to calculate AUC.

# Compute ROC curve
roc_fit1 <- roc(response = mydata$Churn,
                predictor = mydata$EstProb,
                levels = c("No", "Yes"),
                direction = "<")

# Constructs the ROC curve object for model fit1.
# - response: the observed outcome (diabetes status),
# - predictor: predicted probabilities P(Y = 1) from the logistic model,
# - levels: specifies which category is treated as non-event ("neg")
#   and which as event ("pos"),
# - direction: indicates that higher predicted probabilities correspond
#   to a higher likelihood of the event (diabetes).

# Area Under the Curve (AUC)
auc_fit1 <- auc(roc_fit1)
auc_fit1
## Area under the curve: 0.7437
# Computes the area under the ROC curve.
# The AUC summarises the model’s overall discriminatory ability
# independently of any particular classification threshold.

# Plotting ROC curve, including AUC
plot(roc_fit1,
     col = "blue",
     lwd = 2,
     main = paste("ROC Curve for Model fit1 (AUC =",
                  round(auc_fit1, 3), ")"))
# Plots the ROC curve.
# - The y-axis represents sensitivity (true positive rate),
# - The x-axis represents 1 - specificity (false positive rate),
# - The AUC value is displayed in the plot title for interpretation.

# Add reference line for random classification
abline(a = 0, b = 1, lty = 2, col = "gray")

# Adds the 45-degree reference line corresponding to random classification.
# A model with no discriminatory power would lie close to this line.
coords(roc_fit1,
       x = "best",
       best.method = "youden",
       ret = c("threshold", "sensitivity", "specificity", "accuracy"))
##           threshold sensitivity specificity accuracy
## threshold 0.2223175   0.6401571   0.7232199   0.7063